#ifndef BinData2_H
#define BinData2_H

#include <complex>
#include "dbg.h"
#include "Cell.h"

// Some C++ libraries do a fairly complicated thing to multiply two complex numbers
// to make sure the nan and inf's are dealt with correctly.  We don't need that.
// This function does the product as simply as possible.
template <typename T>
inline std::complex<T> SimpleProd(const std::complex<T>& z1, const std::complex<T>& z2)
{
    return std::complex<T>(z1.real()*z2.real()-z1.imag()*z2.imag(),
                           z1.real()*z2.imag()+z1.imag()*z2.real());
}

template <int M>
struct MetricHelper;

template <>
struct MetricHelper<Flat>
{
    static void ProjectShears(
        const Cell<EData,Flat>& c1, const Cell<EData,Flat>& c2,
        const double dsq, std::complex<double>& e1, std::complex<double>& e2)
    {
        // Project given shears to the line connecting them.
        std::complex<double> cr(c2.getData()->pos - c1.getData()->pos);
        std::complex<double> expm2iarg = conj(SimpleProd(cr,cr))/dsq;
        //xxdbg<<"expm2ibeta = "<<expm2iarg<<std::endl;
        e1 = SimpleProd(c1.getData()->we,expm2iarg);
        e2 = SimpleProd(c2.getData()->we,expm2iarg);
    }
};


template <>
struct MetricHelper<Sphere>
{
    static void ProjectShears(
        const Cell<EData,Sphere>& c1, const Cell<EData,Sphere>& c2,
        const double dsq, std::complex<double>& e1, std::complex<double>& e2)
    {
        // For spherical triangles, it's a bit trickier, since the angles aren't equal.
        const Position<Sphere>& p1 = c1.getData()->pos;
        const Position<Sphere>& p2 = c2.getData()->pos;
        //xxdbg<<"p1 = "<<p1<<", p2 = "<<p2<<std::endl;
        // We need the angle at each point between north and the line connecting the 
        // two points.
        //
        // Use the spherical law of cosines:
        //
        // cos(a) = cos(b) cos(c) + sin(b) sin(c) cos(A)
        //
        // In our case:
        //   a = distance from pole to p1 = Pi/2 - dec1
        //   b = distance from pole to p2 = Pi/2 - dec2
        //   c = the great circle distance between the two points.
        //   A = angle between c and north at p2
        //   B = angle between c and north at p1
        //   C = angle between meridians = ra1 - ra2
        // 
        double x1 = p1.getX(); 
        double y1 = p1.getY(); 
        double z1 = p1.getZ(); 
        double x2 = p2.getX(); 
        double y2 = p2.getY(); 
        double z2 = p2.getZ(); 

        // cos(C) = cos(ra1 - ra2) = cos(ra1)cos(ra2) + sin(ra1)sin(ra2)
        //        = (x1/cos(dec1)) (x2/cos(dec2)) + (y1/cos(dec1)) (y2/cos(dec2))
        //        = (x1 x2 + y1 y2) / (cos(dec1) cos(dec2))
        // sin(C) = sin(ra1 - ra2) = sin(ra1)cos(ra2) - cos(ra1)sin(ra2)
        //        = (y1/cos(dec1)) (x2/cos(dec2)) - (x1/cos(dec1)) (y2/cos(dec2))
        //        = (y1 x2 - x1 y2) / (cos(dec1) cos(dec2))

        // cos(A) = (sin(dec1) - sin(dec2) cos(c)) / (cos(dec2) sin(c))
        //
        // The division is fairly unstable if cos(dec2) or sin(c) is small.
        // And sin(c) is often small, so we want to manipulate this a bit.
        // cos(c) = sin(dec1) sin(dec2) + cos(dec1) cos(dec2) cos(C)
        // cos(A) = (sin(dec1) cos(dec2)^2 - sin(dec2) cos(dec1) cos(dec2) cos(C)) /
        //                (cos(dec2) sin(c))
        //        = (sin(dec1) cos(dec2) - sin(dec2) cos(dec1) cos(C)) / sin(c)
        //        = (sin(dec1) cos(dec2)^2 - sin(dec2) (x1 x2 + y1 y2)) / (cos(dec2) sin(c))
        //        = (z1 (1-z2^2) - z2 (x1 x2 + y1 y2)) / (cos(dec2) sin(c))
        // sin(A) / sin(a) = sin(C) / sin(c)
        // sin(A) = cos(dec1) sin(C) / sin(c)
        //        = (y1 x2 - x1 y2) / (cos(dec2) sin(c))
        //
        // We ignore the denominator for now, and then figure out what it is from
        // sin(A)^2 + cos(A)^2 = 1
        double temp = x1*x2+y1*y2;
        double cosA = z1*(1.-z2*z2) - z2*temp;  // These are unnormalized.
        double sinA = y1*x2 - x1*y2;
        //xxdbg<<"A = "<<atan2(sinA,cosA)*180./M_PI<<std::endl;
        double cosAsq = cosA*cosA;
        double sinAsq = sinA*sinA;
        double normAsq = cosAsq + sinAsq;
        double cos2A = (cosAsq - sinAsq) / normAsq; // These are now correct.
        double sin2A = 2.*sinA*cosA / normAsq;
        //xxdbg<<"2A = "<<atan2(sin2A,cos2A)*180./M_PI<<std::endl;

        // cos(B) = (sin(dec2) - sin(dec1) cos(c)) / (cos(dec1) sin(c))
        //        = (sin(dec2) cos(dec1)^2 - sin(dec1) (x1 x2 + y1 y2)) / (cos(dec1) sin(c))
        //        = (z2 (1-z1^2) - z1 (x1 x2 + y1 y2)) / (cos(dec1) sin(c))
        // sin(B) / sin(b) = sin(C) / sin(c)
        // sin(B) = cos(dec2) sin(C) / sin(c)
        //        = (y1 x2 - x1 y2) / (cos(dec1) sin(c))
        double cosB = z2*(1.-z1*z1) - z1*temp;
        //xxdbg<<"B = "<<atan2(sinA,cosB)*180./M_PI<<std::endl;
        // sinB = sinA in the unnormalized state, so just use that.
        double cosBsq = cosB*cosB;
        double normBsq = cosBsq + sinAsq;
        double cos2B = (cosBsq - sinAsq) / normBsq;
        double sin2B = 2.*sinA*cosB / normBsq;
        //xxdbg<<"2B = "<<atan2(sin2B,cos2B)*180./M_PI<<std::endl;

        // In fact, A,B are not really the angles by which we want to rotate the shears.
        // We really want to rotae by the angle between due _east_ and c, not _north_.
        //
        // exp(2ibeta)  = exp(2i (Pi/2 - B) )
        //              = exp(iPi) * exp(-2iB)
        //              = - exp(-2iB)
        //
        // exp(2ialpha) = exp(2i (A - Pi/2) )
        //              = exp(-iPi) * exp(2iA)
        //              = - exp(2iA)
        //
        // We also want the conjugate of these values.

        std::complex<double> expm2ibeta(-cos2B,-sin2B); 
        std::complex<double> expm2ialpha(-cos2A,sin2A);

        //xxdbg<<"expm2ibeta = "<<expm2ibeta<<std::endl;
        //xxdbg<<"expm2ialpha = "<<expm2ialpha<<std::endl;

        e1 = SimpleProd(c1.getData()->we , expm2ibeta);
        e2 = SimpleProd(c2.getData()->we , expm2ialpha);
    }
};

template <class T>
static bool IsNan(T x) { return (x != x) || !(x*x >= 0); }
template <class T>
static bool IsNan(std::complex<T> z) { return IsNan(real(z)) || IsNan(imag(z)); }

template <int DC1, int DC2> struct BinData2;

template <>
struct BinData2<NData,NData>
{
    BinData2() : npair(0.) {}

    double npair;
};

template <>
struct BinData2<EData,EData> 
{
    BinData2() : xiplus(0.), ximinus(0.), varxi(0.), meanlogr(0.), weight(0.), npair(0.) {}

    template <int M>
    void directProcess11(const Cell<EData,M>& c1, const Cell<EData,M>& c2,
                         const double dsq, const double logr)
    {
        xdbg<<"DirectProcess11: p1 = "<<c1.getData()->pos<<", p2 = "<<c2.getData()->pos<<std::endl;
        //xxdbg<<"dsq = "<<dsq<<std::endl;
        //xxdbg<<"logr = "<<logr<<std::endl;
        //xxdbg<<"raw e1,e2 = "<<c1.getData()->we<<", "<<c2.getData()->we<<std::endl;
        std::complex<double> e1, e2;
        MetricHelper<M>::ProjectShears(c1,c2,dsq,e1,e2);
        //xxdbg<<"e1,e2 = "<<e1<<','<<e2<<std::endl;

        // The complex products e1 e2 and e1 e2* share most of the calculations,
        // so faster to do this manually.
        double e1re2r = e1.real() * e2.real();
        double e1re2i = e1.real() * e2.imag();
        double e1ie2r = e1.imag() * e2.real();
        double e1ie2i = e1.imag() * e2.imag();

        const std::complex<double> ee ( e1re2r - e1ie2i , e1ie2r + e1re2i );
        const std::complex<double> eec( e1re2r + e1ie2i , e1ie2r - e1re2i );
        //xxdbg<<"ee,eec = "<<ee<<" "<<eec<<std::endl;

        const double ww = c1.getData()->w * c2.getData()->w;
        const double npairs = c1.getData()->n * c2.getData()->n;

        xiplus += eec;
        ximinus += ee;
        //xxdbg<<"xi => "<<xiplus<<" "<<ximinus<<std::endl;
        weight += ww;
        meanlogr += ww*logr;
        npair += npairs;

#if 0
        if (IsNan(xiplus) || IsNan(ximinus)) {
            dbg<<"Found nan\n";
            dbgout = &std::cout;
            directProcess11(c1,c2,dsq,logr);
            myerror("Found nan");
        }
#endif
    }

    void finalize(double vare)
    {
        xdbg<<"Finalize: weight = "<<weight<<"  vare = "<<vare<<std::endl;
        xdbg<<"   "<<xiplus<<"  "<<ximinus<<"  "<<meanlogr<<"  "<<varxi<<std::endl;
        if (weight == 0.) {
            xiplus = ximinus = meanlogr = varxi = 0.;
        } else {
            xiplus /= weight;
            ximinus /= weight;
            meanlogr /= weight;
            varxi = vare*vare/npair;
        }
        xdbg<<"=> "<<xiplus<<"  "<<ximinus<<"  "<<meanlogr<<"  "<<varxi<<std::endl;
    }

    std::complex<double> xiplus,ximinus;
    double varxi, meanlogr, weight, npair;
};

template <>
struct BinData2<NData,EData> 
{
    BinData2() : meangammat(0.), vargammat(0.), meanlogr(0.), weight(0.), npair(0.) {}

    std::complex<double> meangammat;
    double vargammat, meanlogr, weight, npair;
};


#endif
